%!TEX TS-program = Sweave2

\documentclass[fleqn, 12pt]{article}
\usepackage{geometry}                 % See geometry.pdf to learn the layout options. There are lots.
\geometry{letterpaper}                  % ... or a4paper or a5paper or ... 
%\geometry{landscape}                % Activate for for rotated page geometry
\usepackage[parfill]{parskip}         % Activate to begin paragraphs with an empty line rather than an indent
\usepackage{graphicx}
\usepackage{amsmath, amsthm, amssymb}
\usepackage{epstopdf}
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}
\usepackage{dcolumn}
\usepackage{multicol}
\usepackage{color,soul}
\usepackage{dcolumn}
\usepackage{rotating}

%Packages and preamble borrowed in part from Bill MacMillan:
\usepackage{natbib}
\usepackage{fancyhdr}
\usepackage{sectsty}
\usepackage{setspace}
\usepackage{fullpage}
\usepackage{lscape}
\usepackage{url}
\pagestyle{fancyplain}
\setlength\textheight{9in}
\setlength\voffset{-24pt}
\setlength\headsep{20pt}
\headheight 15pt
%\allsectionsfont{\normalsize\itshape}					%This command sets the section headers to a smaller font.
\newcommand\independent{\protect\mathpalette{\protect\independenT}{\perp}}
\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu{#1#2}}}

%To create hyperrefs
\usepackage{xcolor}
\usepackage[hyperfootnotes=false]{hyperref}
\definecolor{linkcolour}{HTML}{0000cd}	%dark blue color
%\definecolor{shade}{HTML}{D4D7FE}	%light blue shade
%\definecolor{text1}{HTML}{2b2b2b}	%text is almost black
%\definecolor{headings}{HTML}{701112} 	%dark red
\hypersetup{	colorlinks,breaklinks,
			urlcolor=linkcolour, 
			linkcolor=linkcolour}

\usepackage[T1]{fontenc}
\usepackage{mathptmx}  % This is for Times Roman font
%\usepackage{mathpazo}  % This is for Palatino
\usepackage[noae]{Sweave}
%\usepackage{caption}  % This is to personalize captions...

\usepackage[section]{placeins}





%To generate the title page
\title{Put title of document here}
\author{Joshua Gubler}
\date{}                                           % Activate to display a given date or no date

\begin{document}
\lhead{\em Cross-cutting Cleavages and Civil War}
\chead{}
\rhead{Gubler and Selway}

\begin{titlepage}

\begin{center}
\ \\

\vfill
{\LARGE \bf Horizontal Inequality, Cross-cutting Cleavages, and Civil War} \ \\
\ \\
\ \\
{\Large Replication File}
\ \\
\ \\
\ \\
{\large Joshua Gubler\\Assistant Professor\\Department of Political Science\\Brigham Young University}\footnote{Contact: \href{mailto:jgub@byu.edu}{jgub@byu.edu}.}\\
\ \\
{\large Joel Selway\\Assistant Professor\\Department of Political Science\\Brigham Young University}\footnote{Contact: \href{mailto:joel_selway@byu.edu}{joel$\_$selway@byu.edu}.}\\
\ \\
\today


\end{center}
\ \\
\begin{abstract}

This document contains the R-code we used to generate the results in the paper.  In an effort to maximize transparency, it is compiled in $\LaTeX$ and Sweave.  As such, those interested in replication and already running Sweave can simply download the data file as well as the source file for this document, change the file path designated in the line of code to load the data, and then typeset/generate the .pdf with the results.  Those not running Sweave can cut and paste the code in this document into R to replicate the results.  For those running other stats programs, we have saved the data file in .csv and .dta formats as well.

\end{abstract}
\vfill
\end{titlepage}

%To generate the bulk of the paper...


\section{To load the data and necessary R packages for the paper}

<<echo=true,results=hide>>=
library(multilevel)
library(xtable)
library(apsrtable)
library(ggplot2)
library(lmtest)
library(fields)
library(sandwich)
library(Zelig)
##
load("/Users/jrg27/Dropbox/Selway_Projects/JCR_Oct2011_FinalDocs/jcr_27Oct2010.rdata")  
@

\textit{**Note:} To load these data for replication, you will need to change the directory path in the last line of the preceding code to reflect the location of the data file on your computer.


%\section{Code to generate the dataset we will use for our analyses and which we will make public at jcr}

%This is the code we used to generate the final dataset (jcr_27Oct2010.df) from Selway's master dataset:

<<echo=true,results=hide>>=
#jcr_joel.df <- read.dta("/Users/jrg27/Dropbox/Selway_Projects/JCR_RR/jcr_26Oct2010.dta")  ##This is the complete data from Joel's previous work...
###
#jcr_27Oct2010.df <- as.data.frame(jcr_joel.df$cname)
#jcr_27Oct2010.df$CName <- jcr_joel.df$country
#jcr_27Oct2010.df$Country <- jcr_joel.df$CCode2
#jcr_27Oct2010.df$Year <- jcr_joel.df$Year
#jcr_27Oct2010.df$EIC <- jcr_joel.df$ELInc
#jcr_27Oct2010.df$EGC <- jcr_joel.df$ELGeo
#jcr_27Oct2010.df$ERC <- jcr_joel.df$CC
#jcr_27Oct2010.df$ComboCC <- jcr_joel.df$combo3
#jcr_27Oct2010.df$Income <- jcr_joel.df$lgdpenl1
#jcr_27Oct2010.df$Oil <- jcr_joel.df$Oil
#jcr_27Oct2010.df$Population <- jcr_joel.df$lpopl1
#jcr_27Oct2010.df$Mountains <- jcr_joel.df$lmtnest
#jcr_27Oct2010.df$Instability <- jcr_joel.df$instab
#jcr_27Oct2010.df$EthFrac <- jcr_joel.df$EorL
#jcr_27Oct2010.df$RelFrac <- jcr_joel.df$RFrag
#jcr_27Oct2010.df$PITF <- jcr_joel.df$ERWONSET
#jcr_27Oct2010.df$PRIO <- jcr_joel.df$waronset
#jcr_27Oct2010.df$MEPV <- jcr_joel.df$twMEPV
#jcr_27Oct2010.df$PRIOpeaceyrs <- jcr_joel.df$splinevar2
#jcr_27Oct2010.df$PRIOspline1 <- jcr_joel.df$PRIOspline1
#jcr_27Oct2010.df$PRIOspline2 <- jcr_joel.df$PRIOspline2
#jcr_27Oct2010.df$PRIOspline3 <- jcr_joel.df$PRIOspline3
#jcr_27Oct2010.df$PITFpeaceyrs <- jcr_joel.df$PITFpeaceyrs
#jcr_27Oct2010.df$PITFspline1 <- jcr_joel.df$PITFspline1
#jcr_27Oct2010.df$PITFspline2 <- jcr_joel.df$PITFspline2
#jcr_27Oct2010.df$PITFspline3 <- jcr_joel.df$PITFspline3
#jcr_27Oct2010.df$MEPVpeaceyrs <- jcr_joel.df$MEPVpeaceyrs
#jcr_27Oct2010.df$MEPVspline1 <- jcr_joel.df$MEPV_spline1
#jcr_27Oct2010.df$MEPVspline2 <- jcr_joel.df$MEPV_spline2
#jcr_27Oct2010.df$MEPVspline3 <- jcr_joel.df$MEPV_spline3
#jcr_27Oct2010.df <- jcr_27Oct2010.df[,c(-1)]
#rm(jcr_joel.df)
#write.dta(jcr_27Oct2010.df, "jcr_27Oct2010.dta")
#write.csv(jcr_27Oct2010.df, "jcr_27Oct2010.csv")
#save(jcr_27Oct2010.df, file = "jcr_27Oct2010.rdata")
@






\section{To generate the data for Table 1}


\subsection*{Starting with the first row (PITF), going from left to right}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
count(CName[ComboCC>=0 & PITF>=0])
stats(PITF[ComboCC>=0])[2,]
stats(PITF[ComboCC>=0])[3,]
count(CName[ComboCC>=0 & PITF==1])
table(CName[ComboCC>=0 & PITF==1], PITF[ComboCC>=0 & PITF==1])
stats(EthFrac[ComboCC>=0 & PITF>=0])[2,]
stats(RelFrac[ComboCC>=0 & PITF>=0])[2,]
@

\subsection*{Second row (PRIO), going from left to right}

<<echo=true,results=hide>>=
count(CName[ComboCC>=0 & PRIO>=0])
stats(PRIO[ComboCC>=0])[2,]
stats(PRIO[ComboCC>=0])[3,]
count(CName[ComboCC>=0 & PRIO==1])
table(CName[ComboCC>=0 & PRIO==1], PRIO[ComboCC>=0 & PRIO==1])
stats(EthFrac[ComboCC>=0 & PRIO>=0])[2,]
stats(RelFrac[ComboCC>=0 & PRIO>=0])[2,]
@

\subsection*{Third row (MEPV), going from left to right}

<<echo=true,results=hide>>=
count(CName[ComboCC>=0 & MEPV>=0])
stats(MEPV[ComboCC>=0])[2,]
stats(MEPV[ComboCC>=0])[3,]
count(CName[ComboCC>=0 & MEPV==1])
table(CName[ComboCC>=0 & MEPV==1], MEPV[ComboCC>=0 & MEPV==1])
stats(EthFrac[ComboCC>=0 & MEPV>=0])[2,]
stats(RelFrac[ComboCC>=0 & MEPV>=0])[2,]
@

\subsection*{Starting with the fourth row (PITF), going from left to right}

<<echo=true,results=hide>>=
count(CName[PITF>=0])
stats(PITF)[2,]
stats(PITF)[3,]
count(CName[PITF==1])
table(CName[PITF==1], PITF[PITF==1])
stats(EthFrac[PITF>=0])[2,]
stats(RelFrac[PITF>=0])[2,]
@

\subsection*{Fifth row (PRIO), going from left to right}

<<echo=true,results=hide>>=
count(CName[PRIO>=0])
stats(PRIO)[2,]
stats(PRIO)[3,]
count(CName[PRIO==1])
table(CName[PRIO==1], PRIO[PRIO==1])
stats(EthFrac[PRIO>=0])[2,]
stats(RelFrac[PRIO>=0])[2,]
@

\subsection*{Sixth row (MEPV), going from left to right}

<<echo=true,results=hide>>=
count(CName[MEPV>=0])
stats(MEPV)[2,]
stats(MEPV)[3,]
count(CName[MEPV==1])
table(CName[MEPV==1], MEPV[MEPV==1])
stats(EthFrac[MEPV>=0])[2,]
stats(RelFrac[MEPV>=0])[2,]
detach(jcr_27Oct2010.df)
@

%%The table itself: 

\begin{sidewaystable}[htdp]
\begin{center}
\begin{tabular}{|c c c c c c c c|}
\hline\hline
& & & Authors’ Dataset & & & &\\
\hline		
& Total \# of Countries & Mean & Stdev & \# Countries  & \# Civil Wars & Mean & Mean\\
& & & & Experiencing Civil War & & EthFrac & RelFrac \\
PITF & 80 & 0.02& 0.13 & 35 & 59 & 0.42 & 0.46\\
UCDP/PRIO & 80 & 0.01 & 0.10 & 19 & 31 & 0.42 & 0.47\\
MEPV & 80 & 0.01 & 0.10 & 22 & 31 & 0.42 & 0.46\\
\hline
& & & Full Dataset & & & &\\
\hline		
PITF & 159 & 0.02 & 0.13 & 64 & 110 & 0.36 & 0.44\\
UCDP/PRIO & 157 & 0.01 & 0.11 & 49 & 72 & 0.37 & 0.44\\
MEPV & 156 & 0.01 & 0.11 & 55 & 81 & 0.36 & 0.44\\
\hline\hline
\end{tabular}
\end{center}
\label{tab:datasets}
\caption{Summary statistics of countries in our dataset compared to full datasets}
\end{sidewaystable}








\section{To generate the data for Table 2}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
summary.df <- as.data.frame(stats(EIC)[c(2,3,4,8),1])
summary.df$EIC <- stats(EIC)[c(2,3,4,8),1]
summary.df$ERC <- stats(ERC)[c(2,3,4,8),1]
summary.df$EGC <- stats(EGC)[c(2,3,4,8),1]
summary.df$ComboCC <- stats(ComboCC)[c(2,3,4,8),1]
summary.df$Income <- stats(Income)[c(2,3,4,8),1]
summary.df$Oil <- stats(Oil)[c(2,3,4,8),1]
summary.df$Population <- stats(Population)[c(2,3,4,8),1]
summary.df$Mountains <- stats(Mountains)[c(2,3,4,8),1]
summary.df$Instability <- stats(Instability)[c(2,3,4,8),1]
detach(jcr_27Oct2010.df)
##
summary.tab <- t(summary.df)
summary.tab <- summary.tab[c(-1),]
rm(summary.df)
@


%%The table itself:

<<echo=FALSE,results=tex>>= 
print(xtable(summary.tab, caption = "Summary Statistics", label = "tab:sumstats", table.placement = "htbp", caption.placement = "bottom"))
@








\section{To generate the data for the 2x2x2 table in the paper}


\subsection*{To calculate the necessary stats}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
stats(ERC[PITF>=0])
stats(EIC[PITF>=0])
stats(EGC[PITF>=0])
detach(jcr_27Oct2010.df)
@

\subsubsection*{Code to generate the countries, cell by cell, for the TOP row of the table, left to right}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
table(CName[ERC>.75 & EIC>.87 & EGC>.69], PITF[ERC>.75 & EIC>.87 & EGC>.69])
table(CName[ERC<=.75 & EIC>.87 & EGC>.69], PITF[ERC<=.75 & EIC>.87 & EGC>.69])
table(CName[ERC>.75 & EIC<=.87 & EGC>.69], PITF[ERC>.75 & EIC<=.87 & EGC>.69])
table(CName[ERC<=.75 & EIC<=.87 & EGC>.69], PITF[ERC<=.75 & EIC<=.87 & EGC>.69])
detach(jcr_27Oct2010.df)
@

\subsubsection*{Code to generate the countries, cell by cell, for the BOTTOM row of the table, left to right}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
table(CName[ERC>.75 & EIC>.87 & EGC<=.69], PITF[ERC>.75 & EIC>.87 & EGC<=.69])
table(CName[ERC<=.75 & EIC>.87 & EGC<=.69], PITF[ERC<=.75 & EIC>.87 & EGC<=.69])
table(CName[ERC>.75 & EIC<=.87 & EGC<=.69], PITF[ERC>.75 & EIC<=.87 & EGC<=.69])
table(CName[ERC<=.75 & EIC<=.87 & EGC<=.69], PITF[ERC<=.75 & EIC<=.87 & EGC<=.69])
detach(jcr_27Oct2010.df)
@








\section{To generate the regression results in Table 4}


\subsubsection*{Logit estimations, moving from left-to right in the table (using peace spells and splines):}

\textit{**Note:} We created the splines variables using Tucker's ``btscs'' program in Stata.

<<echo=true,results=hide>>=
Model1s <- glm(PITF ~ EIC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model2s <- glm(PITF ~ EGC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model3s <- glm(PITF ~ ERC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model4s <- glm(PITF ~ ComboCC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model5s <- glm(PITF ~ ComboCC + EthFrac + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model6s <- glm(PITF ~ ComboCC + RelFrac + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model7s <- glm(PITF ~ ComboCC + EthFrac + RelFrac + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model8s <- glm(PRIO ~ ComboCC + PRIOpeaceyrs + PRIOspline1 + PRIOspline2 + PRIOspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model9s <- glm(MEPV ~ ComboCC + MEPVpeaceyrs + MEPVspline1 + MEPVspline2 + MEPVspline3 + Income + Oil + Population + Mountains + Instability, family = binomial(link = "logit"), data = jcr_27Oct2010.df)
@



\subsubsection*{Clustered Standard Errors R-code}

To correct the standard errors to account for panel variation, we wrote the following function (adapted from code originally written by Mahmood Arai) which performs the same function as Stata's ``cluster()'' command.  An updated version of this function as well as code for other helpful R functions can be found at Gubler's academic website.

<<echo=true,results=hide>>=
clse <- function(dat,fm, cluster){
 require(sandwich)
 require(lmtest)
 not <- attr(fm$model,"na.action")
	if( ! is.null(not)){ 
  	 cluster <- cluster[-not]
   	 dat <- dat[-not,]
	}
 with(dat,{
 M <- length(unique(cluster))
 N <- length(cluster)
 K <- fm$rank
 dfc <- (M/(M-1))*((N-1)/(N-K))
 uj <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
 vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
 coeftest(fm, vcovCL)
 }
 )
}
@



\subsubsection*{Estimation results with cluster-corrected standard errors (peace spells and splines):}

<<echo=true,results=verbatim>>=
attach(jcr_27Oct2010.df)
clse(jcr_27Oct2010.df,Model1s,Country)
clse(jcr_27Oct2010.df,Model2s,Country)
clse(jcr_27Oct2010.df,Model3s,Country)
clse(jcr_27Oct2010.df,Model4s,Country)
clse(jcr_27Oct2010.df,Model5s,Country)
clse(jcr_27Oct2010.df,Model6s,Country)
clse(jcr_27Oct2010.df,Model7s,Country)
clse(jcr_27Oct2010.df,Model8s,Country)
clse(jcr_27Oct2010.df,Model9s,Country)
detach(jcr_27Oct2010.df)
@





\section{To generate the risk ratios in Table 5}

To calculate the risk ratios, we use the Zelig program created by Imai et. al.  We first re-estimate the models in Table 4 (using GEE, assuming independence within clusters) and then use Zelig's MLE simulation algorithm to run 1000 simulated predictions of Y given the values of X (min/max) we provide.   Here is the code for this procedure:

\subsubsection*{To calculate RR for EIC using Model 1:}

<<echo=true,results=hide>>=
Model1z <- zelig(formula = PITF ~ EIC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, model = "logit.gee", data = jcr_27Oct2010.df, robust = TRUE, id = "Country") 
Model1EIC.high <- setx(Model1z, EIC = 1)
Model1EIC.low <- setx(Model1z, EIC = .538)
Model1.EIC <- sim(Model1z, x = Model1EIC.high, x1 = Model1EIC.low)
summary(Model1.EIC)
@

\subsubsection*{To calculate RR for EGC using Model 2:}

<<echo=true,results=hide>>=
Model2z <- zelig(formula = PITF ~ EGC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, model = "logit.gee", data = jcr_27Oct2010.df, robust = TRUE, id = "Country") 
Model2EGC.high <- setx(Model2z, EGC = 1)
Model2EGC.low <- setx(Model2z, EGC = .538)
Model2.EGC <- sim(Model2z, x = Model2EGC.high, x1 = Model2EGC.low)
summary(Model2.EGC)
@

\subsubsection*{To calculate RR for ERC using Model 3:}

<<echo=true,results=hide>>=
Model3z <- zelig(formula = PITF ~ ERC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, model = "logit.gee", data = jcr_27Oct2010.df, robust = TRUE, id = "Country") 
Model3ERC.high <- setx(Model3z, ERC = 1)
Model3ERC.low <- setx(Model3z, ERC = .538)
Model3.ERC <- sim(Model3z, x = Model3ERC.high, x1 = Model3ERC.low)
summary(Model3.ERC)
@

\subsubsection*{To calculate RR for ComboCC using Model 4:}

<<echo=true,results=hide>>=
Model4z <- zelig(formula = PITF ~ ComboCC + PITFpeaceyrs + PITFspline1 + PITFspline2 + PITFspline3 + Income + Oil + Population + Mountains + Instability, model = "logit.gee", data = jcr_27Oct2010.df, robust = TRUE, id = "Country")
Model4CC.high <- setx(Model4z, ComboCC = 1)
Model4CC.low <- setx(Model4z, ComboCC = .368)
Model4.CC <- sim(Model4z, x = Model4CC.high, x1 = Model4CC.low)
summary(Model4.CC)
@

\subsubsection*{To calculate RR for Income using Model 4:}

<<echo=true,results=hide>>=
Model4Income.high <- setx(Model4z, Income = 11.11)
Model4Income.low <- setx(Model4z, Income = 3.87)
Model4.Income <- sim(Model4z, x = Model4Income.high, x1 = Model4Income.low)
summary(Model4.Income)
@

\subsubsection*{To calculate RR for Oil using Model 4:}

<<echo=true,results=hide>>=
Model4Oil.high <- setx(Model4z, Oil = 1)
Model4Oil.low <- setx(Model4z, Oil = 0)
Model4.Oil <- sim(Model4z, x = Model4Oil.high, x1 = Model4Oil.low)
summary(Model4.Oil)
@

\subsubsection*{To calculate RR for Population using Model 4:}

<<echo=true,results=hide>>=
Model4Population.high <- setx(Model4z, Population = 14.03)
Model4Population.low <- setx(Model4z, Population = 5.40)
Model4.Population <- sim(Model4z, x = Model4Population.high, x1 = Model4Population.low)
summary(Model4.Population)
@

\subsubsection*{To calculate RR for Mountains using Model 4:}

<<echo=true,results=hide>>=
Model4Mountains.high <- setx(Model4z, Mountains = 4.56)
Model4Mountains.low <- setx(Model4z, Mountains = 0)
Model4.Mountains <- sim(Model4z, x = Model4Mountains.high, x1 = Model4Mountains.low)
summary(Model4.Mountains)
@

\subsubsection*{To calculate RR for Instability using Model 4:}

<<echo=true,results=hide>>=
Model4Instability.high <- setx(Model4z, Instability = 1)
Model4Instability.low <- setx(Model4z, Instability = 0)
Model4.Instability <- sim(Model4z, x = Model4Instability.high, x1 = Model4Instability.low)
summary(Model4.Instability)
@







\section{Time dummies code and code to generate the plots presented in the online supplementary appendix}


\subsection*{Table 4 Estimation results with cluster-corrected standard errors (time dummies):}

\subsubsection*{Code to create the time-dummies for PITF}

<<echo=true,results=hide>>=
peacespellPITF <- sort(unique(jcr_27Oct2010.df$PITFpeaceyrs))
timedummiesPITF <- matrix(NA, nrow = nrow(jcr_27Oct2010.df), ncol = length(peacespellPITF))
#
for (j in 1:length(peacespellPITF)) { 
  timedummiesPITF[,j] <- as.integer(jcr_27Oct2010.df$PITFpeaceyrs == peacespellPITF[j])
}
#
pyPITF <- c(0:54)
pyPITF <- paste("psPITF", pyPITF, sep="")
#
dimnames(timedummiesPITF) <- list(1:nrow(timedummiesPITF),pyPITF)
jcr_27Oct2010.df <- cbind(jcr_27Oct2010.df, timedummiesPITF)
@


\subsubsection*{Code to create the time-dummies for PRIO}

<<echo=true,results=hide>>=
peacespellPRIO <- sort(unique(jcr_27Oct2010.df$PRIOpeaceyrs))
timedummiesPRIO <- matrix(NA, nrow = nrow(jcr_27Oct2010.df), ncol = length(peacespellPRIO))
#
for (j in 1:length(peacespellPRIO)) { 
  timedummiesPRIO[,j] <- as.integer(jcr_27Oct2010.df$PRIOpeaceyrs == peacespellPRIO[j])
}
#
pyPRIO <- c(0:49)
pyPRIO <- paste("psPRIO", pyPRIO, sep="")
#
dimnames(timedummiesPRIO) <- list(1:nrow(timedummiesPRIO),pyPRIO)
jcr_27Oct2010.df <- cbind(jcr_27Oct2010.df, timedummiesPRIO)
@


\subsubsection*{Code to create the time-dummies for MEPV}

<<echo=true,results=hide>>=
peacespellMEPV <- sort(unique(jcr_27Oct2010.df$MEPVpeaceyrs))
timedummiesMEPV <- matrix(NA, nrow = nrow(jcr_27Oct2010.df), ncol = length(peacespellMEPV))
#
for (j in 1:length(peacespellMEPV)) { 
  timedummiesMEPV[,j] <- as.integer(jcr_27Oct2010.df$MEPVpeaceyrs == peacespellMEPV[j])
}
#
pyMEPV <- c(0:54)
pyMEPV <- paste("psMEPV", pyMEPV, sep="")
#
dimnames(timedummiesMEPV) <- list(1:nrow(timedummiesMEPV),pyMEPV)
jcr_27Oct2010.df <- cbind(jcr_27Oct2010.df, timedummiesMEPV)
@


\subsubsection*{Logit estimations, moving from left-to right in the table (using time dummies):}

<<echo=true,results=hide>>=
Model1 <- glm(PITF ~ EIC + Income + Oil + Population + Mountains + Instability + psPITF1 + psPITF2 + psPITF3 + psPITF4 + psPITF5 + psPITF6 + psPITF7 + psPITF8 + psPITF9 + psPITF10 + psPITF11 + psPITF12 + psPITF13 + psPITF14 + psPITF15 + psPITF16 + psPITF17 + psPITF18 + psPITF19 + psPITF20 + psPITF21 + psPITF22 + psPITF23 + psPITF24 + psPITF25 + psPITF26 + psPITF27 + psPITF28 + psPITF29 + psPITF30 + psPITF31 + psPITF32 + psPITF33 + psPITF34 + psPITF35 + psPITF36 + psPITF37 + psPITF38 + psPITF39 + psPITF40 + psPITF41 + psPITF42 + psPITF43 + psPITF44 + psPITF45 + psPITF46 + psPITF47 + psPITF48 + psPITF49 + psPITF50, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model2 <- glm(PITF ~ EGC + Income + Oil + Population + Mountains + Instability + psPITF1 + psPITF2 + psPITF3 + psPITF4 + psPITF5 + psPITF6 + psPITF7 + psPITF8 + psPITF9 + psPITF10 + psPITF11 + psPITF12 + psPITF13 + psPITF14 + psPITF15 + psPITF16 + psPITF17 + psPITF18 + psPITF19 + psPITF20 + psPITF21 + psPITF22 + psPITF23 + psPITF24 + psPITF25 + psPITF26 + psPITF27 + psPITF28 + psPITF29 + psPITF30 + psPITF31 + psPITF32 + psPITF33 + psPITF34 + psPITF35 + psPITF36 + psPITF37 + psPITF38 + psPITF39 + psPITF40 + psPITF41 + psPITF42 + psPITF43 + psPITF44 + psPITF45 + psPITF46 + psPITF47 + psPITF48 + psPITF49 + psPITF50, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model3 <- glm(PITF ~ ERC + Income + Oil + Population + Mountains + Instability + psPITF1 + psPITF2 + psPITF3 + psPITF4 + psPITF5 + psPITF6 + psPITF7 + psPITF8 + psPITF9 + psPITF10 + psPITF11 + psPITF12 + psPITF13 + psPITF14 + psPITF15 + psPITF16 + psPITF17 + psPITF18 + psPITF19 + psPITF20 + psPITF21 + psPITF22 + psPITF23 + psPITF24 + psPITF25 + psPITF26 + psPITF27 + psPITF28 + psPITF29 + psPITF30 + psPITF31 + psPITF32 + psPITF33 + psPITF34 + psPITF35 + psPITF36 + psPITF37 + psPITF38 + psPITF39 + psPITF40 + psPITF41 + psPITF42 + psPITF43 + psPITF44 + psPITF45 + psPITF46 + psPITF47 + psPITF48 + psPITF49 + psPITF50, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model4 <- glm(PITF ~ ComboCC + Income + Oil + Population + Mountains + Instability + psPITF1 + psPITF2 + psPITF3 + psPITF4 + psPITF5 + psPITF6 + psPITF7 + psPITF8 + psPITF9 + psPITF10 + psPITF11 + psPITF12 + psPITF13 + psPITF14 + psPITF15 + psPITF16 + psPITF17 + psPITF18 + psPITF19 + psPITF20 + psPITF21 + psPITF22 + psPITF23 + psPITF24 + psPITF25 + psPITF26 + psPITF27 + psPITF28 + psPITF29 + psPITF30 + psPITF31 + psPITF32 + psPITF33 + psPITF34 + psPITF35 + psPITF36 + psPITF37 + psPITF38 + psPITF39 + psPITF40 + psPITF41 + psPITF42 + psPITF43 + psPITF44 + psPITF45 + psPITF46 + psPITF47 + psPITF48 + psPITF49 + psPITF50, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model5 <- glm(PITF ~ ComboCC + EthFrac + Income + Oil + Population + Mountains + Instability + psPITF1 + psPITF2 + psPITF3 + psPITF4 + psPITF5 + psPITF6 + psPITF7 + psPITF8 + psPITF9 + psPITF10 + psPITF11 + psPITF12 + psPITF13 + psPITF14 + psPITF15 + psPITF16 + psPITF17 + psPITF18 + psPITF19 + psPITF20 + psPITF21 + psPITF22 + psPITF23 + psPITF24 + psPITF25 + psPITF26 + psPITF27 + psPITF28 + psPITF29 + psPITF30 + psPITF31 + psPITF32 + psPITF33 + psPITF34 + psPITF35 + psPITF36 + psPITF37 + psPITF38 + psPITF39 + psPITF40 + psPITF41 + psPITF42 + psPITF43 + psPITF44 + psPITF45 + psPITF46 + psPITF47 + psPITF48 + psPITF49 + psPITF50, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model6 <- glm(PITF ~ ComboCC + RelFrac + Income + Oil + Population + Mountains + Instability + psPITF1 + psPITF2 + psPITF3 + psPITF4 + psPITF5 + psPITF6 + psPITF7 + psPITF8 + psPITF9 + psPITF10 + psPITF11 + psPITF12 + psPITF13 + psPITF14 + psPITF15 + psPITF16 + psPITF17 + psPITF18 + psPITF19 + psPITF20 + psPITF21 + psPITF22 + psPITF23 + psPITF24 + psPITF25 + psPITF26 + psPITF27 + psPITF28 + psPITF29 + psPITF30 + psPITF31 + psPITF32 + psPITF33 + psPITF34 + psPITF35 + psPITF36 + psPITF37 + psPITF38 + psPITF39 + psPITF40 + psPITF41 + psPITF42 + psPITF43 + psPITF44 + psPITF45 + psPITF46 + psPITF47 + psPITF48 + psPITF49 + psPITF50, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model7 <- glm(PITF ~ ComboCC + EthFrac + RelFrac + Income + Oil + Population + Mountains + Instability + psPITF1 + psPITF2 + psPITF3 + psPITF4 + psPITF5 + psPITF6 + psPITF7 + psPITF8 + psPITF9 + psPITF10 + psPITF11 + psPITF12 + psPITF13 + psPITF14 + psPITF15 + psPITF16 + psPITF17 + psPITF18 + psPITF19 + psPITF20 + psPITF21 + psPITF22 + psPITF23 + psPITF24 + psPITF25 + psPITF26 + psPITF27 + psPITF28 + psPITF29 + psPITF30 + psPITF31 + psPITF32 + psPITF33 + psPITF34 + psPITF35 + psPITF36 + psPITF37 + psPITF38 + psPITF39 + psPITF40 + psPITF41 + psPITF42 + psPITF43 + psPITF44 + psPITF45 + psPITF46 + psPITF47 + psPITF48 + psPITF49 + psPITF50, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model8 <- glm(PRIO ~ ComboCC + Income + Oil + Population + Mountains + Instability + psPRIO1 + psPRIO2 + psPRIO3 + psPRIO4 + psPRIO5 + psPRIO6 + psPRIO7 + psPRIO8 + psPRIO9 + psPRIO10 + psPRIO11 + psPRIO12 + psPRIO13 + psPRIO14 + psPRIO15 + psPRIO16 + psPRIO17 + psPRIO18 + psPRIO19 + psPRIO20 + psPRIO21 + psPRIO22 + psPRIO23 + psPRIO24 + psPRIO25 + psPRIO26 + psPRIO27 + psPRIO28 + psPRIO29 + psPRIO30 + psPRIO31 + psPRIO32 + psPRIO33 + psPRIO34 + psPRIO35 + psPRIO36 + psPRIO37 + psPRIO38 + psPRIO39 + psPRIO40 + psPRIO41 + psPRIO42 + psPRIO43 + psPRIO44 + psPRIO45 + psPRIO46 + psPRIO47 + psPRIO48 + psPRIO49, family = binomial(link = "logit"), data = jcr_27Oct2010.df)

Model9 <- glm(MEPV ~ ComboCC + Income + Oil + Population + Mountains + Instability + psMEPV1 + psMEPV2 + psMEPV3 + psMEPV4 + psMEPV5 + psMEPV6 + psMEPV7 + psMEPV8 + psMEPV9 + psMEPV10 + psMEPV11 + psMEPV12 + psMEPV13 + psMEPV14 + psMEPV15 + psMEPV16 + psMEPV17 + psMEPV18 + psMEPV19 + psMEPV20 + psMEPV21 + psMEPV22 + psMEPV23 + psMEPV24 + psMEPV25 + psMEPV26 + psMEPV27 + psMEPV28 + psMEPV29 + psMEPV30 + psMEPV31 + psMEPV32 + psMEPV33 + psMEPV34 + psMEPV35 + psMEPV36 + psMEPV37 + psMEPV38 + psMEPV39 + psMEPV40 + psMEPV41 + psMEPV42 + psMEPV43 + psMEPV44 + psMEPV45 + psMEPV46 + psMEPV47 + psMEPV48 + psMEPV49 + psMEPV50 + psMEPV51 + psMEPV52 + psMEPV53 + psMEPV54, family = binomial(link = "logit"), data = jcr_27Oct2010.df)
@


<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
clse(jcr_27Oct2010.df,Model1,Country)
clse(jcr_27Oct2010.df,Model2,Country)
clse(jcr_27Oct2010.df,Model3,Country)
clse(jcr_27Oct2010.df,Model4,Country)
clse(jcr_27Oct2010.df,Model5,Country)
clse(jcr_27Oct2010.df,Model6,Country)
clse(jcr_27Oct2010.df,Model7,Country)
clse(jcr_27Oct2010.df,Model8,Country)
clse(jcr_27Oct2010.df,Model9,Country)
detach(jcr_27Oct2010.df)
@



\subsection*{To generate the scatter plots at the start of the supplementary appendix}

\subsubsection*{Code for the Ethnic Fractionalization and ERC plot:}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
plot(EthFrac,ERC,type="n")
text(EthFrac,ERC,labels=CName)
detach(jcr_27Oct2010.df)
@

\subsubsection*{Code for the Ethnic Fractionalization and EGC plot:}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
plot(EthFrac,EGC,type="n")
text(EthFrac,EGC,labels=CName)
detach(jcr_27Oct2010.df)
@

\subsubsection*{Code for the Ethnic Fractionalization and EIC plot:}

<<echo=true,results=hide>>=
attach(jcr_27Oct2010.df)
plot(EthFrac,EIC,type="n")
text(EthFrac,EIC,labels=CName)
detach(jcr_27Oct2010.df)
@

\ \\
\ \\

\textit{***A concluding note:} We present the raw results for the robustness checks and PCA analysis (both conducted in Stata) in the online supplementary appendix itself.  As such, we do not duplicate the code and results here.  Please contact the authors with any additional questions.





%%To generate the references pages
%\pagebreak
%\bibliographystyle{apsr_fs}
%\bibliography{/Volumes/jgub/Zotero}
%\pagebreak



\end{document} 